Statistik och dataanalys I

F8: Transformationer och multipel linjär regression

Valentin Zulj

Vad har vi gjort hittills?

Vi har introducerat enkel linjär regression, och beskrivit hur vi kan skatta, tolka, och använda en linjär regressionsmodell
Vi har tittat på olika sätt att utvärdera en modell, och kontrollera modellantaganden (m.h.a. t.ex. \(R^2\) och residualanalys)
Vi har även sagt att linjär regression lämpar sig för linjära samband (bild till höger), men inte för icke-linjära samband (bild till vänster)

Vad vi ska göra nu?

Om vi vill analysera ett samband som inte är linjärt kan vi försöka göra det linjärt med hjälp av en transformation
Vi har tittat på log-transformationer tidigare, men kommer nu bredda diskussionen om transformationer lite

Det kan även vara så att vi vill studera sambanden mellan en responsvariabel \(y\) och flera förklarande variabler, \(x_1, \ldots x_p\)
I sådana fall behöver vi använda multipel linjär regression, som vi också kommer ta upp idag

Transformationer

Transformation till linjäritet

Vi återgår till det icke-linjära sambandet (vänster bild), och försöker transformera det till linjäritet
Med transformationerna \(x \rightarrow \log(x)\) och \(y \rightarrow \sqrt{y}\) får vi fram bilden till höger

Sambandet mellan de transformerade variablerna är någorlunda linjärt, och enkel linjär regression passar mycket bättre nu!

Regression med transformerade variabler

Låt oss säga att
- \(x\) är antalet veckor som en patient haft ett läkemedel för sänkt vilopuls
- \(y\) står för samma patients uppmätta vilopuls (slag/minut)
I figurerna får vi då följande axlar

Regression med transformerade variabler

När vi använder de transformerade variablerna blir modellen inte

\[\widehat{\text{vilopuls}} = b_0 + b_1 \cdot \text{behandlingstid}\]

När vi använder de transformerade variablerna får vi i stället

\[\widehat{\sqrt{\text{vilopuls}}} = b_0 + b_1 \cdot \log(\text{behandlingstid)}\]

Vi kan hitta värden för koefficienterna \(b_0\) och \(b_1\), t.ex. med lm()-funktionen i R, och få ut att

\[\widehat{\sqrt{\text{vilopuls}}} = 8.394 - 0.314 \cdot \log(\text{behandlingstid)}\]

Transformera tillbaka till originalskala

Vi kan nu hitta skattningen \(\widehat{\sqrt{\text{vilopuls}}}\), för något värde på \(\log(\text{behandlingstid})\)
Men: vi är troligtvis inte så intresserade av just \(\sqrt{\text{vilopuls}}\) (vad betyder det ens?), utan vårt mål är att skatta just \(\text{vilopuls}\)
Vi är alltså intresserade av vår ursprungliga variabel

För att estimera våra ursprungliga variabler med hjälp av modellen måste vi
1. Transformera det värde av \(x\) som vi är intresserade av
2. Skatta det transformerade värdet av \(y\)-variabeln med vår modell
3. Transformera tillbaka skattningen till den ursprungliga skalan

Transformera tillbaka till originalskala

För att visa ett exempel använder vi vår modell från tidigare, där \[ \widehat{\sqrt{\text{vilopuls}}} = 8.394 - 0.314 \cdot \log(\text{behandlingstid)} \]

Vi vill skatta vilopulsen för en patient som har tagit läkemedlet i fyra veckor
Med andra ord har vi en patient med \(\text{behandlingstid} = 4\), eller \[ \log(\text{behandlingstid}) = \log(4) = 1.386 \]

Regressionsmodellen ger oss nu skattningen

\[ \widehat{\sqrt{\text{vilopuls}}} = 8.394 - 0.314 \cdot 1.386 = 7.959 \]

Transformera tillbaka till originalskala

Vi vet nu att \(\widehat{\sqrt{\text{vilopuls}}} = 7.959\), och det som återstår är att transformera tillbaka till ursprungsskalan, alltså \(\widehat{\sqrt{\text{vilopuls}}}\) till \(\widehat{\text{vilopuls}}\)
För att bli kvitt kvadratroten kvadrerar vi skattningen, och får \[ \widehat{\text{vilopuls}} = \left(\widehat{\sqrt{\text{vilopuls}}}\right)^2 = 7.959^2 = 63.35 \]
Vi får så dra slutsatsen att vilopulsen för en person som behandlats under fyra veckor skattas till ungefär 63 slag per minut

Transformera tillbaka till originalskala

Tabellen nedan visar hur vi transformera tillbaka till originalskala efter några vanliga transformationer

Transformation_av_responsvariabeln	Transformera_tillbaka
\(y \rightarrow y^2\)	\(\hat{y}=\sqrt{\hat{y}^2}\)
\(y \rightarrow y\)	\(\hat{y}=\hat{y}\)
\(y \rightarrow \sqrt{y}\)	\(\hat{y}=\left(\widehat{\sqrt{y}}\right)^2\)
\(y \rightarrow \log(y)\)	\(\hat{y}=e^{\widehat{\log(y)}}\)

Jämförelse med modell utan transformation

Om vi skattar vilopulsen för många olika behandlingstider kan vi dra en linje mellan våra mellan våra skattningar (höger bild)
Linjen fångar mönstret i observationerna mycket bättre än linjen till vänster, som vi får när vi inte transformerar
Med hjälp av transformationer har vi så lyckats fånga ett icke-linjärt samband med hjälp av enkel linjär regression!

Tukeys cirkel – att välja transformation

Val av transformation

I vårt exempel transformerade vi \(y\) till \(\sqrt{y}\) och \(x\) till \(\log(x)\), utan närmare eftertanke – men varför just dessa transformationer?

Till viss del handlar det om att pröva sig fram, och testa olika transformationer tills någon ger ett relativt linjärt samband
Vi behöver dock inte pröva i blindo, utan specifika transformationer passar ofta specifika typer av samband

Tukeys cirkel hjälper oss att identifiera transformationer som kan vara vettinga, baserat på mönstret i vårt spridningsdiagram

Val av transformation

Om vår graf ser ut som på den här bilden kan vi pröva
- Att transformera \(y\) till \(\sqrt{y}\) eller till \(\log(y)\)
- Att transformera \(x\) till \(\sqrt{x}\) eller till \(\log(x)\)

Val av transformation

Om vår graf ser ut som på den här bilden kan vi pröva
- Att transformera \(y\) till \(\sqrt{y}\) eller till \(\log(y)\)
- Att transformera \(x\) till \(x^2\)

Val av transformation

Om vår graf ser ut som på den här bilden kan vi pröva
- Att transformera \(y\) till \(y^2\)
- Att transformera \(x\) till \(\sqrt{x}\) eller till \(\log(x)\)

Transformationer - att välja transformationer

Om vår graf ser ut som på den här bilden kan vi pröva
- Att transformera \(y\) till \(y^2\)
- Att transformera \(x\) till \(x^2\)

Val av transformation

Sammantaget bildar de fyra graferna en cirkel

Tukeys cirkel

Tukeys cirkel beskriver i kompakt form transformeringar som kan prövas
Varje fjärdedel av cirkeln representerar ett visst icke-linjärt mönster
Pilarna uppåt och nedåt indikerar vilka transformationer att pröva för ett visst mönster

Exempel med Tukeys cirkel

Figurerna visar några exempel på transformationer som vi hade kunnat pröva för vårt data om ett läkemedels effekt på vilopulsen
Den sista transformationen ger den “mest” räta linjen, och därför valde vi den

Multipel linjär regression

Hittills har vi pratat om enkel linjär regression, där modellen innehåller en förklarande variabel \[ \hat{y} = b_0 + b_1 x \]

Vi kan dock tänka oss situationer där flera variabler variabler är relevanta
- Om \(y = \text{månadslön}\), är t.ex. både utbildningsnivå och arbetslivserfarenhet rimliga förklarande variabler
- Om \(y = \text{skörd från en jordgubbsplanta}\), kan både temperatur och näringsnivå i jorden vara rimliga \(x\)

En regressionsmodell med flera förklarande variabel kallas för en multipel linjär regressionsmodell

Multipel linjär regression

En multipel linjär regressionsmodell med två förklaringsvariabler skrivs som

\[ \hat{y} = b_0 + b_1 x_1 + b_2 x_2 \]

I detta fall är \(x_1\) och \(x_2\) två olika förklaringsvariabler

Vi kan t.ex. ha

\[ \widehat{\text{jordgubsskörd}} \; \; = b_0 + b_1 \cdot \text{lufttemperatur} + b_2 \cdot \text{näringsnivå} \]

I allmänhet skriver vi en modell med \(p\) variabler som

\[ \hat y = b_0 + b_1 x_1 + b_2 x_2 + \ldots + b_k x_k \]

Multipel linjär regression

I vårt exempel med bilar och bensinförbrukning har vi hittills förklarat bensinförbruningen med bilens vikt
Vi har dock fler variabler i vårt dataset, bland annat variabeln motorstyka, mätt i antalet hästkrafter
Det visar sig att modellen blir bättre om vi använder motorstyrkan som en ytterligare förklaringsvariabel i modellen
Vår nya modell blir

\[ \widehat{\text{litermil}} = b_0 + b_1 \cdot \text{vikt} + b_2 \cdot \text{hästkrafter} \]

Multipel linjär regression i R

Vi kan använda R för att hitta modellens koefficienter, dvs \(b_0\), \(b_1\) och \(b_2\)
Vi lägger till flera förklarande variabler m.h.a. plustecknet

lmod_cars <- lm(litermil ~ viktton + hp, data = mtcars)
lmod_cars$coefficients

(Intercept)     viktton          hp 
0.149155845 0.598453723 0.001769319

Med lite avrundning får vi

\[ \widehat{\text{litermil}} = 0.149 + 0.598 \cdot \text{vikt} + 0.00177 \cdot \text{hästkrafter} \]

Not: I \(\text{mtcars}\) har jag räknat om \(\text{mpg}\) och \(\text{wt}\) till skalorna liter/mil och kg

Beräkning av koefficienter

För enkel linjär regression visade vi formler för beräkning av \(b_0\) och \(b_1\)
När vi gör multipel linjär regression är blir formlerna betydligt mer komplicerade, och ingår därför inte i denna kurs
Vi följer dock samma princip som innan, och väljer de koefficienter som minimerar

\[ \sum_{i=1}^n e_i^2, \;\;\;\;\;\;\;\; e_i^2 = (y_i - \hat{y})^2 \]

Vi väljer alltså koefficienterna som ger den “bästa” modellen, i meningen att den minimerar summan av de kvadrerade residualerna

Generellt antal förklaringsvariabler

Det finns ingen gräns för hur många förklaringsvariabler vi kan använda i en regressionsmodell
Vi såg tidigare att vi kan formulera den generella modellen

\[ \hat{y} = b_0 + b_1 x_1 + b_2 x_2 + ... + b_k x_k \]

Här har vi \(k\) förklaringsvariabler, där \(k\) kan vara vilket antal som helst ¹

Större modeller är t.ex. vanliga när man vill analysera ett samband mellan två specifika variabler, men ta hänsyn till andra variabler som är relevanta

Val av förklaringsvariabler

Vi hävdade att vi kan få en bättre modell för en bils bensinförbrukning genom att använda båda variablerna vikt och häftkrafter
Finns det något sätt att verifiera att detta stämmer?

Ett enkelt sätt är att jämföra \(R^2\) för de båda modellerna – vi vet att modellen med störst \(R^2\) förklarar störst del av variationen i responsvariabeln
På förra föreläsningen såg vi att \(R^2\) blev 0.7919 när vi använde vikt som enda förklaringsvariabel
Nu ska vi se vad \(R^2\) blir när vi lägger till förklaringsvariabeln hästkrafter

Multipel linjär regression - val av förklaringsvariabler

Vi ser att \(R^2=0.8471\), så modellen med två variabler ser ut att vara bättre

summary(lmod_cars)

Kommentar om \(R^2\)

Vi bör dock ta resultatet ovan med en nypa salt, då \(R^2\) alltid blir större när vi lägger till ytterligare variabler
Detta gäller även om det helt saknas samband mellan den nya förklaringsvariabeln och responsvariabeln
Eftersom \(R^2\) blir större oavsett vilken variabel vi lägger till, betyder större \(R^2\) inte nödvändigtvis att modellen blir bättre

Som alternativ kan vi använda ett justerat \(R^2\) (adjusted R-squared)
\(R^2_{\text{adj}}\) är ett mått som ökar med \(R^2\), men som samtidigt minskar samtidigt med antalet förklaringsvariabler
Om den nya variabeln inte förbättrar modellen tillräckligt mycket, kommer därför \(R^2_{\text{adj}}\) att bli mindre

Definition av \(R^2_{\text{adj}}\)

Justerat \(R^2\) beräknas med formeln

\[R_{\text{adj}}^2 = 1 - \cfrac{(1 - R^2)(n-1)}{n-k-1}\]

Här är \(n\) antalet observationer i våra data, och \(k\) antalet förklaringsvariabler i vår modell
Formeln är inte särskilt intuitiv, men i allmänhet gäller som tidigare att större \(R_{\text{adj}}^2\) indikerar en bättre modell

Val av förklaringsvariabler med \(R^2_{\text{adj}}\)

Vi kan nu använda \(R^2_{\text{adj}}\) för att jämföra våra modeller
- Den ursprungliga modellen använder bara vikt som förklaringsvariabel
- Den nya modellen använder även hästkrafter (hp)
Vi börjar med att räkna ut \(R^2_{\text{adj}}\) för den mindre modellen

lmod_cars1 <- lm(litermil ~ viktton, data=mtcars)
summary_1 <- summary(lmod_cars1)
print(summary_1$adj.r.squared)

[1] 0.7849726

Val av förklaringsvariabler med \(R^2_{\text{adj}}\)

Nu undersöker vi vad vi får för resultat när vi lägger till variabeln hästkrafter

lmod_cars2 <- lm(litermil ~ viktton + hp, data=mtcars)
summary_2 <- summary(lmod_cars2)
print(summary_2$adj.r.squared)

[1] 0.8365429

Vi har alltså

	Bara vikt	Vikt och hp
\(R^2_{\text{adj}}\)	0.79	0.84

Vi ser att \(R^2_{\text{adj}}\) stiger när vi lägger till hästkrafter, vilket indikerar att den nya modellen är bättre

Tolkning av koefficienter

Tolkningen av koefficienterna ändras inte mycket, men det finns en liten skillnad som är viktig att ha med sig
Enkel linjär regression: \(\hat{y}\) ökar med \(b_1\) enheter då \(x\) ökar med en enhet

I multipel linjär regression säger vi att \(\hat{y}\) ökar med \(b_j\) enheter då \(x_j\) ökar med en enhet, givet att värdet på övriga förklaringsvariabler hålls konstant
Den sista delen av tolkningen krävs för att tolkningen ska vara korrekt
Om vi ändrar flera variabler samtidigt kan deras effekter på \(\hat y\) förstärka eller ta ut varandra, och tolkningen av ett enskilt \(b_j\) blir då skev

Om \(b_j\) är negativt får vi en “negativ ökning” i \(\hat{y}\) när \(b_j\) ökar – med andra ord får vi en minskning i \(\hat y\)
Vi har följande koefficienter i vår multipla regressionsmodell

(Intercept)     viktton          hp 
0.149155845 0.598453723 0.001769319

Räkneexempel 1: Vi har två bilar som väger exakt lika mycket, men den ena bilen har 1 hästkraft mer än den andra
- Vår skattning är då att bilen med fler hästkrafter förbrukar 0.00177 liter mer per mil

Räkneexempel 2: Vi har nu två bilar med samma antal hästkrafter, men den ena väger 0.2 ton mindre än den andra.
- Modellen uppskattar då att den lättare bilen förbrukar \(0.2 \cdot 0.598 = 0.1196\) liter per mil mindre än den tyngre bilen

Tolkning av koefficienter

Låt oss titta på ett nytt dataset, som vi kan använda för att undersöka sambandet mellan sparande, inkomst och boendekostnad
Vi är intresserade av att skatta hur mycket individerna i vårt dataset sparar årligen, givet deras inkomst och deras boendekostnader
Vi börjar med att göra en linjär regressionsmodell med boendekostnad som enda förklaringsvariabel

Tolkning av koefficienter

Vi skriver modellen som nedan, och beräknar dess koefficienter i R \[ \widehat{\text{sparande}} = b_0 + b_1 \cdot \text{boendekostnad} \]

lm(sparande ~ boendekostnad)$coefficients

  (Intercept) boendekostnad 
 8206.3659478     0.3750134

Om vi avrundar blir vår regressionsmodell

\[\widehat{\text{sparande}} = 8.206 + 0.375 \cdot \text{boendekostnad}\]

Tolkning av koefficienterna

Vi har modellen \[ \widehat{\text{sparande}} = 8.206 + 0.375 \cdot \text{boendekostnad} \]
För varje ytterligare krona en person lägger på boendet per år uppskattar vår modell att personen sparar ytterligare 0.375 kronor per år

Om person A har en boendekostnad som är 10 000 kronor högre än person B, så uppskattar vi att person A sparar 3 750 kronor mer per år

Tolkning av koefficienter

Vi lägger nu till variabeln inkomst, och får en ny modell

lm(sparande ~ boendekostnad + inkomst)$coefficients

  (Intercept) boendekostnad       inkomst 
 -71.74833135   -0.08899729    0.19553001

\[\widehat{\text{sparande}} = -71.7 {\color{red}{-0.089}} \cdot \text{boendekostnad} + 0.196 \cdot \text{inkomst}\]

I denna modell minskar sparandet med boendekostnaden – för varje extra krona i boendekostnad skattar vi att sparandet minskar med 0.089 kronor

När vi bara hade boendekostnad som förklaringsvariabel uppskattade modellen att sparandet ökade med boendekostnaden
Betyder det att en av modellerna har fel? Vilken av dem i så fall?

Tolkning av koefficienterna

Nej, bara för att modellerna ser ut att motsäga varandra, betyder det inte att någon av dem är fel
De två olika modellerna ger oss olika information

Om vi slumpmässigt väljer två individer, så uppskattar den första modellen att den som har högst boendekostnad också är den som sparar mest

\[\widehat{\text{sparande}} = 8.206 + 0.375 \cdot \text{boendekostnad}\]

Här tar vi inte hänsyn till att individerna kan ha olika inkomst
Kan det vara så att den som bor dyrare också har större inkomst, och därför kan spara mer?

Tolkning av koefficienter

Den andra modellen säger: givet värdet på de övriga variablerna kan vi estimera att den som har högre boendekostnad har ett lägre sparande

\[\widehat{\text{sparande}} = -71.7 - 0.089 \cdot \text{boendekostnad} + 0.196 \cdot \text{inkomst}\]

De övriga variablerna är i denna modellen bara inkomst
Modellen säger alltså: givet en viss inkomst, uppskattar vi att den som har högre boendekostnad har ett lägre sparande

Vi jämför nu individer med samma inkomst

Tolkning av koefficienter

Båda modellerna ger logiska resultat, eftersom de beskriver två olika scenarier

Regression 1:

Vi samlar ihop en grupp helt slumpvis utvalda personer
Om vi jämför de personer som har dyrare bostäder med de personer som har billigare bostäder så har de med dyrare bostäder ofta högre inkomster
Det betyder att de också kan spara mer, trots att de har dyrare boende

Regression 2:

Vi samlar ihop en grupp personer där alla har samma inkomst
Inom den här gruppen jämför vi dem som bor dyrt med dem som bor billigt
De som spenderar mindre på sitt boende har mer pengar över att spara, och därför är deras sparande högre

Tolkning av koefficienterna

I multipel linjär regression säger vi att sambandet mellan \(y\) och ett \(x_j\) är betingat på övriga förklaringsvariabler
Jämför med betingade proportioner från F3, då vi undersökte andelen som drabbades av prostatacancer givet att de sällan eller aldrig åt fisk

Det här understyker än en gång att de samband som vi hittar genom regressionsanalyser inte behöver betyda att vi har kausalitet
Regression 1 visar t.ex. inte att högre boendekostnader leder till större sparande, utan snarare att de som har högre boendekostnader kan spara mer trots de högre boendekostnaderna

Tolkning av koefficienter

Vi tittar på ett exempel till, där vi skattar skattar huspriser (i dollar)
Regressionsmodell 1 har antalet sovrum som enda förklaringsvariabel

\[\widehat{\text{price}} = 338.975 + 40.234 \cdot \text{bedroom}\]

Regressionsmodell 2 inkluderar även boytan i kvadratfot (\(ft^2\)).

\[\widehat{\text{price}} = 308.100 + 135 \cdot \text{living area} - 43.347 \cdot \text{bedroom}\]

Koefficienten för antalet sovrum har gått från positiv till negativ
- Vilka möjliga förklaringar kan man tänka sig?
- Hur ska de två olika modellerna tolkas?

Tolkning av koefficienter

Enligt den första regressionsmodellen kostar hus med många sovrum mer än hus med få sovrum – rimligt då hus med många sovrum bör vara större

Enligt den andra modellen kostar hus med många sovrum mindre än hus med få sovrum, givet husets storlek
Om två hus är lika stora uppskattar vi alltså att det hus som har färre sovrum är dyrare
Det skulle kunna bero på att huset med färre sovrum har större kök, större vardagsrum, osv

Fler än 2 förklaringsvariabler

När vi har fler än två förklaringsvariabler måste vi göra tolkningar givet alla de övriga förklaringsvariablerna

Vi tittar på en modell som förklarar bränsleförbrukning med (1) bilens vikt i ton, (2) antalet hästkrafter (hp) och (3) antalet växlar (gear)

 (Intercept)      viktton           hp         gear 
 0.352356646  0.540075281  0.001964705 -0.039753798

Vi ser att koefficienten för antalet växlar är ungefär \(-0.04\), vad betyder detta?

Om vi har två bilar med samma vikt och samma antal hästkrafter, där den enda skillnaden är att bil 1 har en växel mer än bil 2, så uppskattar vi att bil 2 drar 0.04 liter bensin mindre per mil

Förutsättningar för multipel linjär regression

\(y\) och alla \(x_1, x_2, ..., x_k\) måste vara numeriska variabler
\(y\) måste förhålla sig någorlunda linjärt till vardera förklaringsvariabel
Det bör inte finnas uppenbara outliers, eftersom enstaka outliers kan ha stor påverkan på modellens koefficienter
Residualernas varians bör vara konstant
Residualerna bör vara normalfördelade
Residualernas standardavvikelse räknas ut

\[s_e = \sqrt{\cfrac{\sum e^2}{n-k-1}}\]

Här är \(n\) antalet observationer, och \(k\) antalet förklaringsvariabler

Förutsättningar för multipel linjär regression

Vi kan använda nedan figurer för att se om förutsättningarna är uppfyllda
- Spridningsdiagrammet till vänster visar att residualerna på y-axeln inte tycks förändras med våra predikterade värden på x-axeln
- Histogrammet i mitten visar att residualerna är approximativt normalfördelade
- Normalfördelningsgrafen visar även den att residualerna är approximativt normalfördelade

Förutsättningar för multipel linjär regression

Vi behöver inte ställa orimligt höga krav när vi bedömer om en regressionsmodelle uppfyller antagandena
Så länge en modell ger prediktioner som är bättre än \(\bar{y}\) är modellen användbar i rätt sammanhang
En statistisk modell behöver inte vara perfekt, utan det viktiga är att modellens brister kommuniceras öppet!
Då vet den som mottar informationen att modellens resultat ska tolkas med försiktighet

Credits

Dessa slides skapades av Karl Sigfrid för kursen Statistik och Dataanalys I och har uppdaterats av Oskar Gustafsson och Valentin Zulj